% ************************************************************************
% Simulation for the paper "The Perverse Consequences of Policy Restrictions
% in the Presence of Asymmetric Information"
% ************************************************************************
% SECTION 5 (first half)
% ************************************************************************

path(path,'C:\Users\hortalav\Dropbox\Research\_Submitted\RestrictChoice\RestrictChoice (shared)\Old stuff\Matlab');
fclose('all');
clear
clc

global n step step2 t T_1 T_2 A B p alpha max_i U_1 U_2 t1 t2 T1 T2 lu U_NR_u Delta
global w_average_1 w_average_2 w_marginal w_median_1 w_median_2 abst_1 abst_2
global mu_1 mu_2 sigma_1 sigma_2 
global l_1 l_2 l_marginal_1 l_marginal_2 l_u
global epsilon y t_R

% Grid Size
n=1000; %1000
step = 1/n;
step2 = step;
t = (0.05:step:0.95);
max_i=size(t,2);

%%%%%%%%%%%%%%%%%%%%%% LOOP %%%%%%%%%%%%%%%%%%%%%%%%
step3=36;
Results = zeros (step3,step3);
for xx = 1:step3
    for yy = xx:step3  
% the following code captures the case when sigma_1=0.5 & sigma_2 = 1.3
%for xx = 7:1:7
%    for yy = 32:1:32
% FIRST ITERATION: sigma1 and sigma2 go from 0.1 to 2; mu=8
% SECOND ITERATION: mu1 and mu2 go from 0.5 to 15; sigma=1.5  
iteration_min = 0.1;
iteration_max = 1.5;
it_aux1 = ((xx-1) / (step3-1) * (iteration_max-iteration_min) ) + iteration_min;
it_aux2 = ((yy-1) / (step3-1) * (iteration_max-iteration_min) ) + iteration_min;

% Initialisation of parameters 
A=1;
B=1;
alpha = 0.6;

mu_1 = 8;
mu_2 = 8;

sigma_1 = it_aux1;
sigma_2 = it_aux2;

abst_1 = 0.4;
p=0.5;
epsilon = 1; 
y=0.3;

% Equilibrium labour decision of median voters, equilibrium levels of taxes 
% and redistribution when 'abst_1' level of abstention under f_1
w_average_1 = exp(mu_1+(sigma_1^2)/2);
w_average_2 = exp(mu_2+(sigma_2^2)/2);
w_marginal = logninv(abst_1,mu_1,sigma_1);
abst_2 = logncdf(w_marginal,mu_2,sigma_2);
w_median_1 = logninv((abst_1+1)/2,mu_1,sigma_1);
w_median_2 = logninv((abst_2+1)/2,mu_2,sigma_2);
l_1 = zeros (1,max_i);
l_2 = zeros (1,max_i);
U_1 = zeros (1,max_i);
U_2 = zeros (1,max_i);
for i=1:max_i
    l_1(1,i) = (1-alpha)*(1 + t(1,i)* w_average_1 / w_median_1 *(A*t(1,i)+B));
    l_2(1,i) = (1-alpha)*(1 + t(1,i)* w_average_2 / w_median_2 *(A*t(1,i)+B));
    if l_1(1,i)<0
        l_1(1,i)=0;
    elseif l_1(1,i)>1
        l_1(1,i)=1;
    end
    if l_2(1,i)<0
        l_2(1,i)=0;
    elseif l_2(1,i)>1
        l_2(1,i)=1;
    end
    U_1(1,i) = alpha * log ((1-t(1,i))* w_median_1 * (1-l_1(1,i)) + t(1,i)*(1-t(1,i))*(A*t(1,i)+B)*w_average_1)+ (1-alpha)*log(l_1(1,i));
    U_2(1,i) = alpha * log ((1-t(1,i))* w_median_2 * (1-l_2(1,i)) + t(1,i)*(1-t(1,i))*(A*t(1,i)+B)*w_average_2)+ (1-alpha)*log(l_2(1,i)); 
end
[aux,ind_1]=max(U_1);
[aux,ind_2]=max(U_2);
t1=t(1,ind_1);
t2=t(1,ind_2);
T1=t1*(1-t1)*(A*t1+B)*w_average_1;
T2=t2*(1-t1)*(A*t2+B)*w_average_2;

% t1<=t2 (otwerwise swap variables)
if (t1>t2)
    aux=t1;
    t1=t2;
    t2=aux;
    aux=T1;
    T1=T2;
    T2=aux;
    aux=mu_1;
    mu_1=mu_2;
    mu_2=aux;
    aux=sigma_1;
    sigma_1=sigma_2;
    sigma_2=aux;
    aux=w_average_1;
    w_average_1=w_average_2;
    w_average_2=aux;
    p=1-p;
end


% marginal voter: optimal decision if informed and uninformed; cost of
% being informed that makes him indifferent, utility for marginal voter
l_marginal_1 = (1-alpha)*(1 + t1* w_average_1 / w_marginal *(A*t1+B));
l_marginal_2 = (1-alpha)*(1 + t2* w_average_2 / w_marginal *(A*t2+B));
if l_marginal_1<0
    l_marginal_1=0;
elseif l_marginal_1>1
    l_marginal_1=1;
end
if l_marginal_2<0
    l_marginal_2=0;
elseif l_marginal_2>1
    l_marginal_2=1;
end
if (l_marginal_1 ~= l_marginal_2)
    if (l_marginal_1<l_marginal_2)
            aux3=l_marginal_1;
            aux4=l_marginal_2;
    else
            aux3=l_marginal_2;
            aux4=l_marginal_1;
    end
    l_u = (aux3:step2:aux4);
    max_j = size(l_u,2);
    objective_1 = zeros (1,max_j);
    for j=1:max_j
        aux1= p    * (l_u(1,j)-l_marginal_1)*(l_marginal_2-(1-alpha)*l_u(1,j));
        aux2=(1-p) * (l_u(1,j)-l_marginal_2)*(l_marginal_1-(1-alpha)*l_u(1,j));
        objective_1(1,j)= abs(aux1+aux2);
    end
    [aux,ind_u]=min(objective_1);
    lu=l_u(1,ind_u);
else
    lu=l_marginal_1;
end
U_marg_1= alpha * log ((1-t1)* w_marginal * (1-l_marginal_1) + T1)+ (1-alpha)*log(l_marginal_1);
U_marg_2= alpha * log ((1-t2)* w_marginal * (1-l_marginal_2) + T2)+ (1-alpha)*log(l_marginal_2);
U_marg_u_1= alpha * log ((1-t1)* w_marginal * (1-lu) + T1)+ (1-alpha)*log(lu);
U_marg_u_2= alpha * log ((1-t2)* w_marginal * (1-lu) + T2)+ (1-alpha)*log(lu);
Delta = p*(U_marg_1 - U_marg_u_1) + (1-p)*(U_marg_2 - U_marg_u_2);
U_NR_u = p * U_marg_u_1 + (1-p) * U_marg_u_2;

% save NR variables to write them in text file	
t1_NR=t1;
t2_NR=t2;
T1_NR=T1;
T2_NR=T2;
abst_1_NR=abst_1;
abst_2_NR=abst_2;
w_median_1_NR=w_median_1;
w_median_2_NR=w_median_2;
w_marginal_NR=w_marginal;
l_marginal_1_NR =l_marginal_1;
l_marginal_2_NR =l_marginal_2;
l_u_NR=lu;


clearvars l_1 l_2 
% INTRODUCE RESTRICTION: t1=restriction; compute t2 (marginally increase
% abstention and compute optimal decisions of median under f_2; then check
% whether marginal voter has incentives to acquire information; if not,
% increase further abstention; if yes, finish loop; extra condition: stop
% if new tax rate t2 is below the restriction).
t_R = y * t1 + (1-y) * t2;
t1 = t_R;
T1 =  t1*(1-t1)*(A*t1+B)*w_average_1;
flag = true;
w_marginal_R = w_marginal;
while flag
    w_marginal_R= w_marginal_R + epsilon;
    abst_2 =  logncdf(w_marginal_R,mu_2,sigma_2);
    abst_1 = logncdf(w_marginal_R,mu_1,sigma_1);        %not needed for computation but nice for data
    w_median_1 = logninv((abst_1+1)/2,mu_1,sigma_1); 	%not needed for computation but nice for data
    w_median_2 = logninv((abst_2 + 1)/2,mu_2,sigma_2);
    l_2 = zeros (1,max_i);
    for i=1:max_i
        l_2(1,i) = (1-alpha)*(1 + t(1,i)* w_average_2 / w_median_2 *(A*t(1,i)+B));       
        if l_2(1,i)<0
            l_2(1,i)=0;
        elseif l_2(1,i)>1
            l_2(1,i)=1;
        end      
        U_2(1,i) = alpha * log ((1-t(1,i))* w_median_2 * (1-l_2(1,i)) + t(1,i)*(1-t(1,i))*(A*t(1,i)+B)*w_average_2)+ (1-alpha)*log(l_2(1,i)); 
    end
    [aux,ind_2]=max(U_2);
    t2=t(1,ind_2);
    if (t2<=t_R)
        flag=false;
        t2=t_R;
    end
    T2=t2*(1-t2)*(A*t2+B)*w_average_2;
    l_marginal_1 = (1-alpha)*(1 + t1* w_average_1 / w_marginal_R *(A*t1+B));
    l_marginal_2 = (1-alpha)*(1 + t2* w_average_2 / w_marginal_R *(A*t2+B));
    if l_marginal_1<0
        l_marginal_1=0;
    elseif l_marginal_1>1
        l_marginal_1=1;
    end
    if l_marginal_2<0
        l_marginal_2=0;
    elseif l_marginal_2>1
        l_marginal_2=1;
    end    
    
    
    if (l_marginal_1 ~= l_marginal_2)
        if (l_marginal_1<l_marginal_2)
            aux3=l_marginal_1;
            aux4=l_marginal_2;
        else
            aux3=l_marginal_2;
            aux4=l_marginal_1;
        end
        l_u_R = (aux3:step2:aux4);
        max_j = size(l_u_R,2);
        objective_2 = zeros (1,max_j);
        for j=1:max_j
            aux1= p    * (l_u_R(1,j)-l_marginal_1)*(l_marginal_2-(1-alpha)* l_u_R(1,j));
            aux2=(1-p) * (l_u_R(1,j)-l_marginal_2)*(l_marginal_1-(1-alpha)* l_u_R(1,j));
            objective_2(1,j)= abs(aux1+aux2);
        end
        [aux,ind_u]=min(objective_2);
        lu= l_u_R(1,ind_u);
    else
          lu=l_marginal_1;
    end
    U_marg_1= alpha * log ((1-t1)* w_marginal_R * (1-l_marginal_1) + T1)+ (1-alpha)*log(l_marginal_1);
    U_marg_2= alpha * log ((1-t2)* w_marginal_R * (1-l_marginal_2) + T2)+ (1-alpha)*log(l_marginal_2);
    U_marg_u_1= alpha * log ((1-t1)* w_marginal_R * (1-lu) + T1)+ (1-alpha)*log(lu);
    U_marg_u_2= alpha * log ((1-t2)* w_marginal_R * (1-lu) + T2)+ (1-alpha)*log(lu);
    Delta_R = p*(U_marg_1 - U_marg_u_1) + (1-p)*(U_marg_2 - U_marg_u_2);
%    if (Delta_R >= Delta) || (w_marginal_R > w_marginal + 1000000*epsilon)
    if (Delta_R >= Delta)
        flag=false;
    end
end


% FINAL: compute current welfare of no-restriction-marginal-voter
l_marginal_1_cost = (1-alpha)*(1 + t1* w_average_1 / w_marginal *(A*t1+B));
l_marginal_2_cost = (1-alpha)*(1 + t2* w_average_2 / w_marginal *(A*t2+B));
if l_marginal_1_cost<0
    l_marginal_1_cost=0;
elseif l_marginal_1_cost>1
    l_marginal_1_cost=1;
end
if l_marginal_2_cost<0
    l_marginal_2_cost=0;
elseif l_marginal_2_cost>1
    l_marginal_2_cost=1;
end
if (l_marginal_1_cost ~= l_marginal_2_cost)
    if (l_marginal_1_cost<l_marginal_2_cost)
        aux3=l_marginal_1_cost;
        aux4=l_marginal_2_cost;
    else
        aux3=l_marginal_2_cost;
        aux4=l_marginal_1_cost;
    end
    l_u_R_cost = (aux3:step2:aux4);
    max_j = size(l_u_R_cost,2);
    objective_3 = zeros (1,max_j);
    for j=1:max_j
        aux1= p    * (l_u_R_cost(1,j)-l_marginal_1_cost)*(l_marginal_2_cost-(1-alpha)* l_u_R_cost(1,j));
        aux2=(1-p) * (l_u_R_cost(1,j)-l_marginal_2_cost)*(l_marginal_1_cost-(1-alpha)* l_u_R_cost(1,j));
        objective_3(1,j)= abs(aux1+aux2);
    end
    [aux,ind_u]=min(objective_3);
    lu_cost= l_u_R_cost(1,ind_u);
else
      lu_cost=l_marginal_1_cost;
end
U_marg_u_1= alpha * log ((1-t1)* w_marginal * (1-lu_cost) + T1)+ (1-alpha)*log(lu_cost);
U_marg_u_2= alpha * log ((1-t2)* w_marginal * (1-lu_cost) + T2)+ (1-alpha)*log(lu_cost);
U_R_u = p * U_marg_u_1 + (1-p) * U_marg_u_2;
% Cost_of_Restriction = U_NR_u - U_R_u;  % when positive, restriction has unintended consequences
Cost_of_Restriction = alpha*(p*log( T1_NR)+(1-p)*log( T2_NR)-p*log( T1)-(1-p)*log( T2));  % when positive, restriction has unintended consequences

auxx = [num2str(xx),' , ',num2str(yy),' , ' num2str(Cost_of_Restriction)];
disp(auxx)
Results(xx,yy)=Cost_of_Restriction;
Results(yy,xx)=Cost_of_Restriction; %it is symmetric but we only execute the triangle.
end
end

% plot matrix of results
figure
mesh(Results)
title({'Loss when introducing the restriction';'(as we vary the standard deviation of income distributions)'})
title({'Loss when introducing the restriction' ; ' for citizens that subsist on welfare payments' ;'(as we vary the standard deviation of income distributions)'})

xlabel('sigma_1')
ylabel('sigma_2')

ax = gca;
set(ax,'XTick',[1  floor(step3/3) floor(2*step3/3) step3])
set(ax,'YTick',[1  floor(step3/3) floor(2*step3/3) step3])
tick1=iteration_min;
tick2 = ((floor(step3/3)-1) / (step3-1) * (iteration_max-iteration_min) ) + iteration_min;
tick3 = ((floor(2*step3/3)-1) / (step3-1) * (iteration_max-iteration_min) ) + iteration_min;
tick4=iteration_max;

set(ax,'Ztick',[])
set(ax,'XTickLabel',{tick1, 0.5, 1, tick4})
set(ax,'YTickLabel',{tick1, 0.5, 1, tick4})
